# ---------------------------
# R version: 4.0.2
# Code: 9/18/2020
# Title: 'Police Demilitarization and Violent Crime'
# Author: Kenneth Lowande
# Summary: Replication code for producing all tables and figures in the main text; replication code and additional data for the Supplementary Information will be provided unconditionally, upon request. Please report errors to: lowande@umich.edu
# ---------------------------

cat("\014")
rm(list=ls())
# install.packages(c('lfe','ggplot2','reshape2','grid','gridExtra','ggthemes'))
require(lfe) # v2.8-5
require(ggplot2) # v3.3.2
require(reshape2) # v1.4.4
require(grid) # v4.0.2
require(gridExtra) # v2.3
require(ggthemes) # v4.2.0

setwd('...')

# Figure 1:  Roughly one in five controlled weapons turnover within 4 years

# load 2014-Q4 shipments and differences
load('q4-2014.Rda')
load('differences.Rda')

# list NSNs of items to subset
assault_rifles=c('1005-00-073-9421','1005-00-856-6885','1005-01-128-9936','1005-01-042-9820','1005-00-933-7672','6920-01-513-6665','1005-00-589-1271')
shotguns=c('1005-01-073-2368','1005-00-677-9150','1005-00-921-5483')
light_armor=c('2355-21-913-4649','2355-01-123-1606','2320-01-123-1601','2355-01-481-8580')
mine_resistant=c('2355-01-554-8556','2355-01-553-4634','2355-01-562-6146','2355-01-602-3357','2355-01-623-0967','2355-01-555-0908','2355-01-590-1660','2355-01-562-6152','2355-01-596-1330','2355-01-552-4677','2355-01-570-7996','2355-01-561-0281','2355-01-552-5199','2355-01-529-2246','2355-01-552-5344','2355-01-553-4636','2355-01-621-5506','2355-01-552-5173','2355-01-564-4340','2355-01-558-1053','2355-01-562-6675','2355-01-552-5581')
fast_attack='2320-01-467-0677'
handguns=c('1005-00-726-5655','1005-01-D17-2240')

# tracking attrition of items
track_attrition=function(items,data,change,label){
  qtrs=c("2015-1","2015-2","2015-3","2015-4","2016-1","2016-2","2016-3","2016-4","2017-1","2017-2","2017-3","2017-4","2018-1","2018-2","2018-3","2018-4","2019-1","2019-2","2019-3","2019-4","2020-1")
  start=data[data$NSN%in%items,]
  attrition=change[change$NSN%in%items&change$change.type=='Attrition',]
  result=data.frame(qtr=c('2014-4',qtrs),value=numeric(22),number=numeric(22),perc=numeric(22),item=rep(label,22))
  result[1,2]=sum(start$Quantity)
  result[1,3]=sum(start$Acquisition.Value)
  result[1,4]=1
  temp_inv=start
  for (i in 1:length(qtrs)){
    temp=attrition[attrition$qtr==qtrs[i],]
    temp_inv=temp_inv[!temp_inv$uid%in%temp$uid,]
    result[1+i,2]=sum(temp_inv$Quantity)
    result[1+i,3]=sum(temp_inv$Acquisition.Value)
    result[1+i,4]=sum(temp_inv$Quantity)/result[1,2]
  }
  return(result)
}

# summarize by type
hg=track_attrition(handguns,q14,diffs,'Handguns')
sg=track_attrition(shotguns,q14,diffs,'Shotguns')
ar=track_attrition(assault_rifles,q14,diffs,'Assault Rifles')
mr=track_attrition(mine_resistant,q14,diffs,'MRVs')
AT=rbind(hg,sg,ar,mr)
AT=AT[AT$qtr=='2015-1'|AT$qtr=='2017-1'|AT$qtr=='2019-4',]

# plot time snapshots
ATTR=ggplot(AT,aes(x=qtr,y=perc,group=item,color=item)) + 
  theme_bw() + xlab('') + ylab('% Remaining') + geom_line() + geom_point() + ylim(0,1) +
  scale_colour_manual(values = c("#000000","#7F0000","#989898", "#00274c")) + labs(color = "Dec. 2014 Inventory") +
  scale_x_discrete(breaks=c('2015-1','2017-1','2019-4'),labels=c(2015,2017,2019)) +
  theme(panel.background=element_blank())
# ---------------------------

# Figure  2:  Demilitarization  did  not  lead  to  substantively  significant  increases  or  decreases  in violent crime or officer safety

# load lists of treated, matched agencies, and weights
load('treated-matched.Rda')

# load complete and generate treatment variable
load('main.Rda')

# subset to 6 quarters
dt=dt[dt$qtr!='2014-4' & dt$qtr!='2015-1' & dt$qtr!='2016-4' & dt$qtr!='2017-1' & dt$qtr!='2017-2' & dt$qtr!='2017-3' & dt$qtr!='2017-4' & dt$qtr!='2018-1' & dt$qtr!='2018-2' & dt$qtr!='2018-3' & dt$qtr!='2018-4',]

# subset data to this list and add weights - baseline 
dtm1=dt[dt$ori9%in%mat1$IDS,]
dtm1=merge(dtm1,data.frame(ori9=mat1$IDS,weights=mat1$weights),by='ori9')
vars=c('actual_murder','treated','ori9','qtr')
dtm1=dtm1[complete.cases(dtm1[,vars]),]

# subset data to this list and add weights - restrictive 
dtm2=dt[dt$ori9%in%mat2$IDS,]
dtm2=merge(dtm2,data.frame(ori9=mat2$IDS,weights=mat2$weights),by='ori9')
vars=c('actual_murder','treated','ori9','qtr')
dtm2=dtm2[complete.cases(dtm2[,vars]),]

# subset data to this list and add weights - inclusive 
dtm3=dt[dt$ori9%in%mat3$IDS,]
dtm3=merge(dtm3,data.frame(ori9=mat3$IDS,weights=mat3$weights),by='ori9')
vars=c('actual_murder','treated','ori9','qtr')
dtm3=dtm3[complete.cases(dtm3[,vars]),]

# labels
l=c('Murder','Assault','Rape','Manslaughter','Robbery','Felonious Officer Deaths','Accidental Officer Deaths','Officer Assaults (Injury)','Officer Assaults (No Injury)')

# formulas
f=list(as.formula(log(actual_murder+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(actual_assault_total+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(actual_rape_total+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(actual_manslaughter+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(actual_robbery_total+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(officers_killed_by_felony+1) ~ treated| as.factor(ori9) + as.factor(qtr)),
       as.formula(log(officers_killed_by_accident+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(assaults_with_injury_total+1) ~ treated | as.factor(ori9) + as.factor(qtr)),
       as.formula(log(assaults_no_injury_total+1) ~ treated | as.factor(ori9) + as.factor(qtr)))

# estimation and results storage
results=data.frame(sample=character(length(l)*3),outcome=character(length(l)*3),b=numeric(length(l)*3),se=numeric(length(l)*3),stringsAsFactors = F)

for (i in 1:length(f)){
  m=felm(f[[i]],dtm1,weights=dtm1$weights,na.action=na.omit)
  results[i,'b']=as.numeric(m$coefficients)
  results[i,'se']=as.numeric(m$se)
  results[i,'outcome']=l[i] 
  results[i,'sample']='Baseline'
  
  m=felm(f[[i]],dtm2,weights=dtm2$weights,na.action=na.omit)
  results[i+9,'b']=as.numeric(m$coefficients)
  results[i+9,'se']=as.numeric(m$se)
  results[i+9,'outcome']=l[i] 
  results[i+9,'sample']='Restrictive'
  
  m=felm(f[[i]],dtm3,weights=dtm3$weights,na.action=na.omit)
  results[i+18,'b']=as.numeric(m$coefficients)
  results[i+18,'se']=as.numeric(m$se)
  results[i+18,'outcome']=l[i] 
  results[i+18,'sample']='Inclusive'
}

# plot results
interval1 <- -qnorm((1-0.9)/2) 
interval2 <- -qnorm((1-0.95)/2)
interval3 <- -qnorm((1-0.99)/2)
results$outcome=factor(results$outcome,levels=rev(l))

F_RES=ggplot(results,aes(colour = sample)) + scale_colour_manual(values = c("#989898","#00274c","#7F0000"), name = 'Sample') +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  geom_linerange(aes(x = outcome, ymin = b - se*interval1,ymax = b + se*interval1),
                 lwd = 3/2,position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x = outcome, y = b, ymin = b-se*interval2,ymax = b+se*interval2),
                  lwd = 1, position = position_dodge(width = 1/2),shape = 21) +
  geom_linerange(aes(x = outcome, ymin = b - se*interval3,ymax = b + se*interval3),
                 lwd = 1/2,position = position_dodge(width = 1/2)) +
  coord_flip() + theme_bw() + xlab('') + ylab('')
# ---------------------------

# Figure 3:  EO 13688 achieved near-perfect compliance and resulted in demilitarization among affected  LEAs
# NOTE: Also produces Supplementary Figures 1-2, which appear in the Supplementary Information
load('summarized-inventory-2014-2020.Rda')

x=summaryBy(attr2_val+tran2_val+new2_val~qtr,data=d_,FUN=sum,keep.names=T)
x=melt(x,id.vars='qtr')
levels(x$variable)=c('Attrition','Transfers','New Distributions')

F_TS_A=ggplot(x,aes(x=qtr,y=value/1000000,group=variable,color=variable)) + 
  theme_bw() + xlab('') + ylab('Value (Millions)') + geom_line(size=1.5) +
  scale_colour_manual(values = c("#7F0000","#989898", "#00274c","#000000")) + labs(color = "") +
  scale_x_discrete(breaks=c('2015-1', '2016-1', '2017-1','2018-1','2019-1','2020-1'),labels=c(2015, 2016,2017,2018,2019,2020)) +
  theme(panel.grid.major=element_blank(),panel.border=element_blank(),
        panel.background=element_blank())

x=summaryBy(attr2_quant+tran2_quant+new2_quant~qtr,data=d_,FUN=sum,keep.names=T)
x=melt(x,id.vars='qtr')
levels(x$variable)=c('Attrition','Transfers','New Distributions')

F_TS_B=ggplot(x,aes(x=qtr,y=value/1000,group=variable,color=variable)) + 
  theme_bw() + xlab('') + ylab('Number (Thousands)') + geom_line(size=1.5) +
  scale_colour_manual(values = c("#7F0000","#989898", "#00274c","#000000")) + labs(color = "") +
  scale_x_discrete(breaks=c('2015-1', '2016-1', '2017-1','2018-1','2019-1','2020-1'),labels=c(2015, 2016,2017,2018,2019,2020)) +
  theme(panel.grid.major=element_blank(),panel.border=element_blank(),
        panel.background=element_blank())

x=summaryBy(inv_r2_val+attr_r2_val+tran_r2_val+new_r2_val~qtr,data=d_,FUN=sum,keep.names=T)
x=melt(x,id.vars='qtr')
levels(x$variable)=c('In Circulation','Attrition','Transfers','New Distributions')

F_RE_A=ggplot(x,aes(x=qtr,y=value/1000000,group=variable,color=variable)) + 
  theme_bw() + xlab('') + ylab('Value (Millions)') + geom_line(size=1.5) +
  scale_colour_manual(values = c("#000000","#7F0000","#989898", "#00274c")) + labs(color = "") +
  annotate("rect", xmin = 4, xmax = 5, ymin = 0, ymax = 34,alpha = .45)+
  scale_x_discrete(breaks=c('2015-1','2016-1','2017-1','2018-1','2019-1','2020-1'),labels=c(2015,2016,2017,2018,2019,2020)) +
  theme(panel.grid.major=element_blank(),panel.border=element_blank(),
        panel.background=element_blank())

x=summaryBy(inv_r2_ships+attr_r2_ships+tran_r2_ships+new_r2_ships~qtr,data=d_,FUN=sum,keep.names=T)
x=melt(x,id.vars='qtr')
levels(x$variable)=c('In Circulation','Attrition','Transfers','New Distributions')

F_RE_B=ggplot(x,aes(x=qtr,y=value,group=variable,color=variable)) + 
  theme_bw() + xlab('') + ylab('Shipments') + geom_line(size=1.5) +
  annotate("rect", xmin = 4, xmax = 5, ymin = 0, ymax = 600,alpha = .45)+
  scale_colour_manual(values = c("#000000","#7F0000","#989898", "#00274c")) + labs(color = "") +
  scale_x_discrete(breaks=c('2015-1','2016-1','2017-1','2018-1','2019-1','2020-1'),labels=c(2015,2016,2017,2018,2019,2020)) +
  theme(panel.grid.major=element_blank(),panel.border=element_blank(),
        panel.background=element_blank())

g_legend=function(a.gplot){
  tmp=ggplot_gtable(ggplot_build(a.gplot))
  leg=which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend=tmp$grobs[[leg]]
  return(legend)}
legend=g_legend(F_RE_A)

grid.arrange(arrangeGrob(F_RE_A + theme(legend.position="none"),nrow=1),arrangeGrob(F_RE_B + theme(legend.position="none")),legend,nrow=1,widths=c(2,2,1))
# ---------------------------

# Figure 4:  Local  law  enforcement  agencies  with  recalled  equipment  were  concentrated  in  the South, Midwest, and California
# NOTE: Also produces Supplementary Figure 13, which appears in the Supplementary Information
load('usa-shape.Rda')

# aggregate pre-treatment equipment data from Dec 2014
load('q14-pretreatment.Rda')
d=summaryBy(inv_c_val+inv_r2_val~id,data=dt_pre,FUN=sum,keep.names=T)
d=merge(d,data.frame(id=us$id),by='id',all.y=T)
d=d[!duplicated(d),]
d$inv_c_val[is.na(d$inv_c_val)]=0
d$inv_r2_val[is.na(d$inv_r2_val)]=0

# map of all controlled equipment
MAP1=ggplot() + geom_map(data = d, aes(map_id = id, fill = log(inv_c_val+1)),map = us) +
  geom_polygon(data=us, aes(x=long, y=lat, group = group), size = 0.25, colour='white', fill=NA) +
  expand_limits(x = us$long, y = us$lat) + 
  theme_map() + theme(legend.position = 'bottom') +
  scale_fill_gradient(name = 'Log(Equipment Value): ',low='#e5e9ed',high='#00274c')

# map of recalled equipment
MAP2=ggplot() + geom_map(data = d, aes(map_id = id, fill = log(inv_r2_val+1)),map = us) +
  geom_polygon(data=us, aes(x=long, y=lat, group = group), size = 0.25, colour='white', fill=NA) +
  expand_limits(x = us$long, y = us$lat) + 
  theme_map() + theme(legend.position = 'bottom') +
  scale_fill_gradient(name = 'Log(Recalled Equipment Value): ',low='#f2e5e5',high='#7F0000')

# load pre-processed matches lat-long file
load('matches-latlong.Rda')

# map that locates matched agencies
MAP3=ggplot() + geom_polygon(data=us, aes(x=long, y=lat, group = group), size = 0.25, colour='grey', fill=NA) +
  geom_point(data = mat.plot, aes(x = long, y = lat, size = `LEA Group`, shape = `LEA Group`, color = `LEA Group`, alpha = `LEA Group`)) +
  scale_colour_manual(values = c('#00274c','#7F0000')) + 
  scale_size_manual(values = c(1.8,3)) +
  scale_alpha_manual(values = c(0.4,1)) +
  expand_limits(x = us$long, y = us$lat) + 
  theme_map() + theme(legend.position = c(.07,.05))
# ---------------------------
